home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / msdos / fractal / frasr182 / mpmath_c.c < prev    next >
C/C++ Source or Header  |  1993-04-16  |  12KB  |  465 lines
  1. /* MPMath_c.c (C) 1989, Mark C. Peterson, CompuServe [70441,3353]
  2.      All rights reserved.
  3.  
  4.    Code may be used in any program provided the author is credited
  5.      either during program execution or in the documentation.  Source
  6.      code may be distributed only in combination with public domain or
  7.      shareware source code.  Source code may be modified provided the
  8.      copyright notice and this message is left unchanged and all
  9.      modifications are clearly documented.
  10.  
  11.      I would appreciate a copy of any work which incorporates this code,
  12.      however this is optional.
  13.  
  14.      Mark C. Peterson
  15.      405-C Queen St. Suite #181
  16.      Southington, CT 06489
  17.      (203) 276-9721
  18. */
  19.  
  20. #include <stdlib.h>
  21. #include "mpmath.h"
  22. #include "prototyp.h"
  23.  
  24. #ifndef XFRACT
  25.  
  26. int MPaccuracy = 32;
  27.  
  28. struct MP *MPsub(struct MP x, struct MP y) {
  29.    y.Exp ^= 0x8000;
  30.    return(MPadd(x, y));
  31. }
  32.  
  33. /* added by TW */
  34. struct MP *MPsub086(struct MP x, struct MP y) {
  35.    y.Exp ^= 0x8000;
  36.    return(MPadd086(x, y));
  37. }
  38.  
  39. /* added by TW */
  40. struct MP *MPsub386(struct MP x, struct MP y) {
  41.    y.Exp ^= 0x8000;
  42.    return(MPadd386(x, y));
  43. }
  44.  
  45. struct MP *MPabs(struct MP x) {
  46.    x.Exp &= 0x7fff;
  47.    return(&x);
  48. }
  49.  
  50. struct MPC MPCsqr(struct MPC x) {
  51.    struct MPC z;
  52.  
  53.         z.x = *pMPsub(*pMPmul(x.x, x.x), *pMPmul(x.y, x.y));
  54.         z.y = *pMPmul(x.x, x.y);
  55.         z.y.Exp++;
  56.    return(z);
  57. }
  58.  
  59. struct MP MPCmod(struct MPC x) {
  60.         return(*pMPadd(*pMPmul(x.x, x.x), *pMPmul(x.y, x.y)));
  61. }
  62.  
  63. struct MPC MPCmul(struct MPC x, struct MPC y) {
  64.    struct MPC z;
  65.  
  66.         z.x = *pMPsub(*pMPmul(x.x, y.x), *pMPmul(x.y, y.y));
  67.         z.y = *pMPadd(*pMPmul(x.x, y.y), *pMPmul(x.y, y.x));
  68.    return(z);
  69. }
  70.  
  71. struct MPC MPCdiv(struct MPC x, struct MPC y) {
  72.    struct MP mod;
  73.  
  74.    mod = MPCmod(y);
  75.         y.y.Exp ^= 0x8000;
  76.         y.x = *pMPdiv(y.x, mod);
  77.         y.y = *pMPdiv(y.y, mod);
  78.    return(MPCmul(x, y));
  79. }
  80.  
  81. struct MPC MPCadd(struct MPC x, struct MPC y) {
  82.    struct MPC z;
  83.  
  84.         z.x = *pMPadd(x.x, y.x);
  85.         z.y = *pMPadd(x.y, y.y);
  86.    return(z);
  87. }
  88.  
  89. struct MPC MPCsub(struct MPC x, struct MPC y) {
  90.    struct MPC z;
  91.  
  92.         z.x = *pMPsub(x.x, y.x);
  93.         z.y = *pMPsub(x.y, y.y);
  94.    return(z);
  95. }
  96.  
  97. struct MPC MPCone = { 0x3fff, 0x80000000l, 0, 0l };
  98.  
  99. struct MPC MPCpow(struct MPC x, int exp) {
  100.    struct MPC z;
  101.    struct MPC zz;
  102.  
  103.    if(exp & 1)
  104.       z = x;
  105.    else
  106.       z = MPCone;
  107.    exp >>= 1;
  108.    while(exp) {
  109.                 zz.x = *pMPsub(*pMPmul(x.x, x.x), *pMPmul(x.y, x.y));
  110.                 zz.y = *pMPmul(x.x, x.y);
  111.                 zz.y.Exp++;
  112.       x = zz;
  113.       if(exp & 1) {
  114.                         zz.x = *pMPsub(*pMPmul(z.x, x.x), *pMPmul(z.y, x.y));
  115.                         zz.y = *pMPadd(*pMPmul(z.x, x.y), *pMPmul(z.y, x.x));
  116.          z = zz;
  117.       }
  118.       exp >>= 1;
  119.    }
  120.    return(z);
  121. }
  122.  
  123. int MPCcmp(struct MPC x, struct MPC y) {
  124.    struct MPC z;
  125.  
  126.         if(pMPcmp(x.x, y.x) || pMPcmp(x.y, y.y)) {
  127.                 z.x = MPCmod(x);
  128.                 z.y = MPCmod(y);
  129.                 return(pMPcmp(z.x, z.y));
  130.    }
  131.    else
  132.       return(0);
  133. }
  134.  
  135. _CMPLX MPC2cmplx(struct MPC x) {
  136.    _CMPLX z;
  137.  
  138.         z.x = *pMP2d(x.x);
  139.         z.y = *pMP2d(x.y);
  140.    return(z);
  141. }
  142.  
  143. struct MPC cmplx2MPC(_CMPLX z) {
  144.    struct MPC x;
  145.  
  146.         x.x = *pd2MP(z.x);
  147.         x.y = *pd2MP(z.y);
  148.    return(x);
  149. }
  150.  
  151. /* function pointer versions added by Tim Wegner 12/07/89 */
  152. int        (*ppMPcmp)() = MPcmp086;
  153. int        (*pMPcmp)(struct MP x, struct MP y) = MPcmp086;
  154. struct MP  *(*pMPmul)(struct MP x, struct MP y)= MPmul086;
  155. struct MP  *(*pMPdiv)(struct MP x, struct MP y)= MPdiv086;
  156. struct MP  *(*pMPadd)(struct MP x, struct MP y)= MPadd086;
  157. struct MP  *(*pMPsub)(struct MP x, struct MP y)= MPsub086;
  158. struct MP  *(*pd2MP)(double x)                 = d2MP086 ;
  159. double *(*pMP2d)(struct MP m)                  = MP2d086 ;
  160. struct MP  *(*pfg2MP)(long x, int fg)          = fg2MP086;
  161.  
  162. void setMPfunctions(void) {
  163.    if(cpu == 386)
  164.    {
  165.       pMPmul = MPmul386;
  166.       pMPdiv = MPdiv386;
  167.       pMPadd = MPadd386;
  168.       pMPsub = MPsub386;
  169.       pMPcmp = MPcmp386;
  170.       pd2MP  = d2MP386 ;
  171.       pMP2d  = MP2d386 ;
  172.       pfg2MP = fg2MP386;
  173.    }
  174.    else
  175.    {
  176.       pMPmul = MPmul086;
  177.       pMPdiv = MPdiv086;
  178.       pMPadd = MPadd086;
  179.       pMPsub = MPsub086;
  180.       pMPcmp = MPcmp086;
  181.       pd2MP  = d2MP086 ;
  182.       pMP2d  = MP2d086 ;
  183.       pfg2MP = fg2MP086;
  184.    }
  185. }
  186.  
  187. #endif /* XFRACT */
  188.  
  189. #ifndef sqr
  190. #define sqr(x) ((x)*(x))
  191. #endif
  192.  
  193. extern int debugflag, fpu;
  194.  
  195. _CMPLX ComplexPower(_CMPLX xx, _CMPLX yy) {
  196.    _CMPLX z, cLog, t;
  197.    double e2x, siny, cosy;
  198.  
  199.    FPUcplxlog(&xx, &cLog);
  200.    FPUcplxmul(&cLog, &yy, &t);
  201.  
  202.    if(fpu == 387)
  203.       FPUcplxexp387(&t, &z);
  204.    else {
  205.       if(t.x < -690)
  206.          e2x = 0;
  207.       else
  208.          e2x = exp(t.x);
  209.       FPUsincos(&t.y, &siny, &cosy);
  210.       z.x = e2x * cosy;
  211.       z.y = e2x * siny;
  212.    }
  213.    return(z);
  214. }
  215.  
  216. /***** FRACTINT specific routines and variables *****/
  217.  
  218. #ifndef TESTING_MATH
  219.  
  220. #include <stdlib.h>
  221.  
  222. #include "fractint.h"
  223.  
  224. extern double param[];
  225. extern _CMPLX old, new, init;
  226. extern double threshold, roverd, d1overd, dx0[], dy0[];
  227. extern int periodicitycheck, row, col, debugflag;
  228.  
  229.  
  230. #include <stdlib.h>
  231.  
  232. extern int  xdots, ydots;     /* coordinates of dots on the screen  */
  233. extern int  colors;           /* maximum colors available */
  234. extern int  maxit;
  235.  
  236. BYTE far *LogTable = (BYTE far *)0;
  237. extern int LogFlag;
  238.    /* LogFlag == 1  -- standard log palettes
  239.       LogFlag == -1 -- 'old' log palettes
  240.       LogFlag >  1  -- compress counts < LogFlag into color #1
  241.       LogFlag < -1  -- use quadratic palettes based on square roots && compress
  242.    */
  243.  
  244. void SetupLogTable(void) {
  245.    float l, f, c, m;
  246.    unsigned n, prev, limit, lf;
  247.  
  248.    if (LogFlag > -2) {
  249.       lf = (LogFlag > 1) ? LogFlag : 0;
  250.       if (lf >= maxit)
  251.          lf = maxit - 1;
  252.       Fg2Float((long)(maxit-lf), 0, m);
  253.       fLog14(m, m);
  254.       Fg2Float((long)(colors-(lf?2:1)), 0, c);
  255.       fDiv(m, c, m);
  256.       for (prev = 1; prev <= lf; prev++)
  257.          LogTable[prev] = 1;
  258.       for (n = (lf?2:1); n < colors; n++) {
  259.          Fg2Float((long)n, 0, f);
  260.          fMul16(f, m, f);
  261.          fExp14(f, l);
  262.          limit = Float2Fg(l, 0) + lf;
  263.          if (limit > maxit || n == colors-1)
  264.             limit = maxit;
  265.          while (prev <= limit)
  266.             LogTable[prev++] = n;
  267.       }
  268.    } else {
  269.       if ((lf = 0 - LogFlag) >= maxit)
  270.          lf = maxit - 1;
  271.       Fg2Float((long)(maxit-lf), 0, m);
  272.       fSqrt14(m, m);
  273.       Fg2Float((long)(colors-2), 0, c);
  274.       fDiv(m, c, m);
  275.       for (prev = 1; prev <= lf; prev++)
  276.          LogTable[prev] = 1;
  277.       for (n = 2; n < colors; n++) {
  278.          Fg2Float((long)n, 0, f);
  279.          fMul16(f, m, f);
  280.          fMul16(f, f, l);
  281.          limit = Float2Fg(l, 0) + lf;
  282.          if (limit > maxit || n == colors-1)
  283.             limit = maxit;
  284.          while (prev <= limit)
  285.             LogTable[prev++] = n;
  286.       }
  287.    }
  288.    LogTable[0] = 0;
  289.    if (LogFlag != -1)
  290.       for (n = 1; n < maxit; n++) /* spread top to incl unused colors */
  291.          if (LogTable[n] > LogTable[n-1])
  292.             LogTable[n] = LogTable[n-1]+1;
  293. }
  294.  
  295. long far ExpFloat14(long xx) {
  296.    static float fLogTwo = (float)0.6931472;
  297.    int f;
  298.    long Ans;
  299.  
  300.    f = 23 - (int)RegFloat2Fg(RegDivFloat(xx, *(long*)&fLogTwo), 0);
  301.    Ans = ExpFudged(RegFloat2Fg(xx, 16), f);
  302.    return(RegFg2Float(Ans, (char)f));
  303. }
  304.  
  305. extern _CMPLX tmp;
  306. extern int color, colors;
  307. double TwoPi;
  308. _CMPLX temp, t2, BaseLog;
  309. _CMPLX cdegree = { 3.0, 0.0 },
  310.                croot   = { 1.0, 0.0 };
  311.  
  312. int ComplexNewtonSetup(void) {
  313.    threshold = .001;
  314.    periodicitycheck = 0;
  315.    if(param[0] != 0.0 || param[1] != 0.0 || param[2] != 0.0 ||
  316.       param[3] != 0.0) {
  317.       croot.x = param[2];
  318.       croot.y = param[3];
  319.       cdegree.x = param[0];
  320.       cdegree.y = param[1];
  321.       FPUcplxlog(&croot, &BaseLog);
  322.       TwoPi = asin(1.0) * 4;
  323.    }
  324.    return(1);
  325. }
  326.  
  327. int ComplexNewton(void) {
  328.    _CMPLX cd1;
  329.  
  330.    /* new = ((cdegree-1) * old**cdegree) + croot
  331.             ----------------------------------
  332.                  cdegree * old**(cdegree-1)         */
  333.  
  334.    cd1.x = cdegree.x - 1.0;
  335.    cd1.y = cdegree.y;
  336.  
  337.    temp = ComplexPower(old, cd1);
  338.    FPUcplxmul(&temp, &old, &new);
  339.  
  340.    tmp.x = new.x - croot.x;
  341.    tmp.y = new.y - croot.y;
  342.    if((sqr(tmp.x) + sqr(tmp.y)) < threshold)
  343.       return(1);
  344.  
  345.    FPUcplxmul(&new, &cd1, &tmp);
  346.    tmp.x += croot.x;
  347.    tmp.y += croot.y;
  348.  
  349.    FPUcplxmul(&temp, &cdegree, &t2);
  350.    FPUcplxdiv(&tmp, &t2, &old);
  351.    if(DivideOverflow)
  352.    {
  353.       DivideOverflow = 0;
  354.       return(1);
  355.    }
  356.    new = old;
  357.    return(0);
  358. }
  359.  
  360. int ComplexBasin(void) {
  361.    _CMPLX cd1;
  362.    double mod;
  363.  
  364.    /* new = ((cdegree-1) * old**cdegree) + croot
  365.             ----------------------------------
  366.                  cdegree * old**(cdegree-1)         */
  367.  
  368.    cd1.x = cdegree.x - 1.0;
  369.    cd1.y = cdegree.y;
  370.  
  371.    temp = ComplexPower(old, cd1);
  372.    FPUcplxmul(&temp, &old, &new);
  373.  
  374.    tmp.x = new.x - croot.x;
  375.    tmp.y = new.y - croot.y;
  376.    if((sqr(tmp.x) + sqr(tmp.y)) < threshold) {
  377.       if(fabs(old.y) < .01)
  378.          old.y = 0.0;
  379.       FPUcplxlog(&old, &temp);
  380.       FPUcplxmul(&temp, &cdegree, &tmp);
  381.       mod = tmp.y/TwoPi;
  382.       color = (int)mod;
  383.       if(fabs(mod - color) > 0.5) {
  384.          if(mod < 0.0)
  385.             color--;
  386.          else
  387.             color++;
  388.       }
  389.       color += 2;
  390.       if(color < 0)
  391.          color += 128;
  392.       return(1);
  393.    }
  394.  
  395.    FPUcplxmul(&new, &cd1, &tmp);
  396.    tmp.x += croot.x;
  397.    tmp.y += croot.y;
  398.  
  399.    FPUcplxmul(&temp, &cdegree, &t2);
  400.    FPUcplxdiv(&tmp, &t2, &old);
  401.    if(DivideOverflow)
  402.    {
  403.       DivideOverflow = 0;
  404.       return(1);
  405.    }
  406.    new = old;
  407.    return(0);
  408. }
  409.  
  410. extern int Distribution, Offset, Slope;
  411. extern long con;
  412.  
  413. /*** PB, commented this out, it was unused, actual work is in prompts.c
  414. int Starfield(void) {
  415.    int c;
  416.  
  417.    plasma();
  418.    Distribution = (int)param[1];
  419.    con = (long)(param[2] / 100 * (1L << 16));
  420.    Slope = (int)param[3];
  421.    for(row = 0; row < ydots; row++) {
  422.       for(col = 0; col < xdots; col++) {
  423.          if(check_key())
  424.             return(-1);
  425.          c = getcolor(col, row);
  426.          putcolor(col, row, GausianNumber(c, colors));
  427.       }
  428.    }
  429.    return(0);
  430. }
  431.   ***/
  432.  
  433. /*
  434.  * Generate a gaussian distributed number.
  435.  * The right half of the distribution is folded onto the lower half.
  436.  * That is, the curve slopes up to the peak and then drops to 0.
  437.  * The larger slope is, the smaller the standard deviation.
  438.  * The values vary from 0+offset to range+offset, with the peak
  439.  * at range+offset.
  440.  * To make this more complicated, you only have a
  441.  * 1 in Distribution*(1-Probability/Range*con)+1 chance of getting a
  442.  * Gaussian; otherwise you just get offset.
  443.  */
  444. int GausianNumber(int Probability, int Range) {
  445.    int n, r;
  446.    long Accum = 0, p;
  447.  
  448.    p = divide((long)Probability << 16, (long)Range << 16, 16);
  449.    p = multiply(p, con, 16);
  450.    p = multiply((long)Distribution << 16, p, 16);
  451.    if(!(rand15() % (Distribution - (int)(p >> 16) + 1))) {
  452.       for(n = 0; n < Slope; n++)
  453.          Accum += rand15();
  454.       Accum /= Slope;
  455.       r = (int)(multiply((long)Range << 15, Accum, 15) >> 14);
  456.       r = r - Range;
  457.       if(r < 0)
  458.          r = -r;
  459.       return(Range - r + Offset);
  460.    }
  461.    return(Offset);
  462. }
  463.  
  464. #endif
  465.